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Abstract 


Although the study of X-ray binaries has led to major breakthroughs in high-energy astrophysics, their circumbinary 
environment at scales of —100—10,000 au has not been thoroughly investigated. In this paper, we undertake a novel 
and exploratory study by employing direct and high-contrast imaging techniques on a sample of X-ray binaries, using 
adaptive optics and the vortex coronagraph on Keck/NIRC2. High-contrast imaging opens up the possibility to 
search for exoplanets, brown dwarfs, circumbinary companion stars, and protoplanetary disks in these extreme 
systems. Here we present the first near-infrared high-contrast images of 13 high-mass X-ray binaries located within 
~2-3 kpc. The key results of this campaign involve the discovery of several candidate circumbinary companions 
ranging from substellar (brown dwarf) to stellar masses. By conducting an analysis based on Galactic population 
models, we discriminate sources that are likely background/foreground stars and isolate those that have a high 
probability (26096-9996) of being gravitationally bound to the X-ray binary. This paper seeks to establish a 
preliminary catalog for future analyses of proper motion and subsequent observations. With our preliminary results, 
we calculate the first estimate of the companion frequency and the multiplicity frequency for X-ray binaries: 20.6 
and 1.8 + 0.9, respectively, considering only the sources that are most likely bound to the X-ray binary. In addition to 
extending our comprehension of how brown dwarfs and stars can form and survive in such extreme systems, our 
study opens a new window to our understanding of the formation of X-ray binaries. 


Unified Astronomy Thesaurus concepts: Multiple stars (1081); X-ray binary stars (1811); Near infrared astronomy 
(1093); High mass x-ray binary stars (733); Coronagraphic imaging (313); High contrast techniques (2369); Direct 
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imaging (387); Exoplanet detection methods (489); Substellar companion stars (1648) 


1. Introduction 


X-ray binaries are semidetached binary systems in which a 
compact object (white dwarf (WD), neutron star (NS), or stellar- 
mass black hole (BH)) accretes material from a donor star. These 
systems undergo several extreme physical phenomena, such as 
processes acting predominantly in soft X-rays (e.g., Khargharia 
et al. 2010; Tetarenko et al. 2021) and detectable X-ray 
pulsations (e.g., Lutovinov et al. 2005). 

The variations in physical processes among different X-ray 
binaries are directly linked to the mass of the donor star. Over 
90% of these systems can be classified into two distinct 
categories: high-mass X-ray binaries (HMXBs; Maonor Z 8 Mo) 
and low-mass X-ray binaries (LMXBs; Maonor S 1, 5 Mo; e.g. 
Tauris & van den Heuvel 2006). LMXBs are relatively old 
systems (>10° yr) harboring a K-M spectral type donor star, 
where the process of mass transfer is believed to be triggered by 
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Roche lobe overflow (RLO; e.g., Savonije 1978). RLO is 
triggered either by stellar evolution or by angular momentum 
loss (e.g., Paczyński 1967; Verbunt & Zwaan 1981; Strohmayer 
2002; Justham et al. 2006; Chen & Podsiadlowski 2016; Seto 
2018; Van et al. 2019). The transferred mass then agglomerates 
to form an accretion disk around the compact object, giving rise 
to transient accretion and X-ray emission (e.g., Charles & 
Coe 2006). 

As for HMXBs, they are generally thought to be younger 
systems (<10’ yr) harboring a massive O-B spectral type donor 
star. The transferred and accreted matter is thought to 
predominantly come from the capture of a fraction of the stellar 
winds ejected from the donor star (e.g Mohamed & 
Podsiadlowski 2007; Abate et al. 2013; El Mellah et al. 2019). 
There are two subcategories of HMXBs relevant to this work. 
First, we emphasize Be/X-ray binaries (BeXRBs), wherein the 
donor star is a fast-rotating Be star. In these systems, the X-ray 
emission is mainly triggered by the compact object passing 
through a diffuse and gaseous circumstellar disk surrounding the 
Be star (known as a decretion disk; e.g., Okazaki et al. 2002; 
Martin et al. 2011; Rimulo et al. 2018; Kravtsov et al. 2020). 
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Second, we highlight supergiant fast X-ray transients (SFXTs; 
Negueruela et al. 2006), characterized by the presence of a 
supergiant donor star and by fast transient X-ray flaring activity 
within the system (likely induced by an NS; e.g., Sidoli 2012; 
Ducci et al. 2019). 

X-ray binaries are important touchstone objects for high- 
energy phenomena in astrophysics. They have been widely used 
to study several high-energy astronomical phenomena, including 
accretion physics (e.g., Done et al. 2007; Kara et al. 2019) and 
outflow/jet processes (e.g., Markoff et al. 2001; Fender et al. 
2004; Mooley et al. 2018). However, the immediate surround- 
ings of X-ray binaries, at the scale of ~100—10,000 astronomical 
units (au), have been poorly studied. This paper undertakes a 
pioneering exploration of the circumstellar environments of 
X-ray binaries through the application of adaptive optics (AO) 
and direct/high-contrast imaging techniques. The goal is to 
probe a variety of phenomena, ranging from protoplanetary 
disks to debris disks and fallback disks, and particularly to 
search for wide-orbiting circumbinary companions (CBCs)—be 
they exoplanets, brown dwarfs, or stars. 

Considering the discovery of planetary-mass CBCs orbiting 
both binary systems (e.g., Bakos et al. 2007; Desidera & 
Barbieri 2007; Eriksson et al. 2020) and compact objects (WDs 
or pulsars; e.g., Wolszczan & Frail 1992; Sigurdsson et al. 
2003; Spiewak et al. 2018; Vanderburg et al. 2020; Blackman 
et al. 2021), it is not unfounded to expect CBCs orbiting X-ray 
binaries. A recent study argued that X-ray binaries could host 
planetary systems in close orbits detectable via X-ray eclipses 
(Imara & Di Stefano 2018). In this paper, we explore wider 
orbits (~100-10,000 au), as the increased number of interac- 
tions within the system could lead to the ejection of potential 
CBCs from the direct environment of the X-ray binary (e.g., 
Bonavita et al. 2016). 

In Prasow-Emond et al. (2022), we presented the first set of 
observations from a pilot study aiming to survey all X-ray 
binaries amenable for direct imaging within ~3 kpc. We first 
targeted a y Cassiopeiae—like X-ray binary harboring a Be 
donor star, RX J1744.7—2713, for which we had observations 
from two different bands and two epochs. We unveiled the 
presence of three potential CBCs within this system, exhibiting 
a strong likelihood of being stellar-mass CBCs. Here we 
present the first L’-band high-contrast images of 13 other 
systems and conduct a preliminary statistical analysis derived 
from the results of the first epochs of observations. 

The paper is organized as follows: Section 2 presents the 
sample and how it was constructed. Section 3 presents the near- 
infrared observations and the data reduction and processing. 
Section 5 presents the first high-contrast images of the observed 
X-ray binaries. In Section 6, we analyze the images and explore 
the nature of the detection. Finally, in Section 7, we discuss our 
results and their implications. 


2. The Sample 


Despite the ongoing active search for new X-ray binaries 
both within and beyond our Galaxy (e.g., Gandhi et al. 2022), 
their presence remains relatively scarce. Our Galaxy hosts 
~300 identified X-ray binaries known to date (Liu et al. 2006, 
2007). We drew on this list of X-ray binaries as our initial data 
set; however, not all of these systems are suitable for direct 
imaging with Keck/NIRC2. In order to build a sample of X-ray 
binaries that would yield optimal statistical constraints and 
mitigate potential biases, we used four selection criteria: 
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. Distance. The system must be close enough to resolve the 
direct environment at ~100—10,000 au scales. We chose a 
distance limit of ~3 kpc within our Galaxy, which enables 
the detection of structures and objects located within a 
couple of thousands of au from the X-ray binary. The outer 
limit (10,000 au) corresponds to the approximate limit of 
the Keck/NIRC2 field of view (FOV). 

2. Brightness and Adaptive Optics. The donor star must be 
bright enough (Z< 9-10 mag) for the AO loop to be 
closed. 

3. Age. We targeted young (100 Myr) X-ray binaries to 
favor the detection of substellar CBCs, considering the 
steep decline in planet brightness with time (e.g., 
Burrows et al. 2001). This limited us to bright LMXBs 
and HMXBs with a massive O/B donor star. 

4. Visibility. We selected X-ray binaries visible from the W. 

M. Keck Observatory at the time of our observations 

(Keck Observatory Semesters 2017B and 2020A). 


Applying these criteria narrowed down the initial list from 
~300 X-ray binaries to 19, out of which 14 were observed 
between 2017 and 2020 using Keck/NIRC2 (see Section 3.1 
for more details). Our sample includes both HMXBs and 
LMXBs (e.g., MAXI J1820+070, V404 Cyg, 1A 0620-00); 
however, we have only observed HMXBs to date owing to 
observational constraints. Note that X-ray binary surveys are 
far from complete, as X-ray binaries can be undetectable in 
quiescence (e.g., Bird et al. 2007; Belczynski & Ziolkowski 
2009). New X-ray binaries have also been discovered since we 
did our sampling (Avakyan et al. 2023; Neumann et al. 2023). 
Nonetheless, we assume that our sample is as complete as 
possible and that our results accurately represent the known 
population of X-ray binaries. 

Figure 1 displays the position of the 14 observed targets onto 
the sky (Aitoff projection). It also indicates the distance from the 
observer (inkpc) and, if known, the nature of the system's 
compact object. In Section 4, we present a brief literature review 
for each of the observed X-ray binaries. Table 1 summarizes the 
known relevant physical properties of the systems, namely the 
X-ray binary type, X-ray emission class, donor star spectral type, 
and compact object type. Additional relevant physical properties 
can be found in the Appendix (see Table 5). 


3. Observations and Data Reduction 
3.1. Keck/NIRC2 Observations 


On 2017 September 8, we observed four HMXBs from 
our survey sample using the Keck/NIRC2 vortex corona- 
graph (Mawet et al. 2005; Serabyn et al. 2016) in pupil-tracking 
mode in L’ band (A = 3.776 um, AA = 0.700 um; PI: Mawet), 
and with the narrow camera (plate scale of 9.971 + 0.004 mas 
pixel; Service et al. 2016). On 2018 January 3, we observed 
three additional HMXBs using a similar setting. Due to the 
successful and promising preliminary results of this initial 
campaign, we were awarded three supplementary nights of 
observation on 2020 July 11, 12, and 13 (PI: Fogarty). On the 
first night, we observed one additional target and reobserved two 
targets (RX J1744.7—2713 and ^ Cas) using a similar setting. 
However, due to saturation, we had to downscale the frame size 
of y Cas from 1024 x 1024 pixel? to 512 x 512 pixel’. On 2020 
July 12, we obtained data for three other HMXBs in L’ band, in 
addition to reobserving RX J1744.7—2713 in K, band 
(A= 2.146 um, AA = 0.311 um). Finally, on 2020 July 13, we 
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Figure 1. Position of the 14 observed X-ray binaries on the sky with equatorial coordinates and Aitoff projection, color-coded with distance from the observer in kpc. 
The nature of the compact object is illustrated with different markers: squares for WDs, crosses for NSs, asterisks for BHs, and circles for cases where the nature is 
either unavailable or uncertain. The red dotted line shows the approximate coordinates of the Galactic plane. 


obtained data for three other HMXBs, totaling observations for 
14 out of the 19 X-ray binaries in the sample. 

During the observations, we used the Quadrant Analysis of 
Coronagraphic Images for Tip-tilt Sensing (Huby et al. 2017) 
to make tip-tilt adjustments to maintain precise centering of the 
target on the vortex focal plane mask. The observations were 
AO assisted using the Shack-Hartmann wave front sensor 
(which performs wave front sensing in R band) in 2017, in 
2018, and on the last night of 2020. For the first two nights of 
our 2020 observations, we opted for the Pyramid wave front 
sensor instead. It performs wave front sensing in H 
(Wizinowich et al. 2000; Bond et al. 2018), which is more 
advantageous for the redder targets in our sample. 

A summary of the observing log is presented in Table 2. 


3.2. Data Reduction 


Similarly to Prasow-Émond et al. (2022), we performed data 
reduction using the Vortex Image Processing (VIP) and NIRC2 
Preprocessing packages (Gomez Gonzalez et al. 2017). To 
obtain a preprocessed data cube, we proceeded as follows: (1) 
flat-fielding of the frames, (2) bad pixel masking using the dark 
frames, (3) determination of the vortex center for each frame 
followed by cropping the science cube around the mean center, 
(4) removal of sky contribution via a principal component 
analysis (PCA) based technique, and (5) image registration to 
align the quasi-static speckle pattern across frames. 

After acquiring the preprocessed data cube, we applied a 
PCA-based angular differential imaging (ADI; Marois et al. 
2006) algorithm to obtain high-contrast images. Subsequently, 
we generated several images using two algorithms in VIP 
(annular PCA and full-frame PCA) for a broad range of 
principal components (from 1 to 50). This was done to ensure 
consistency in the detection of sources within the images (i.e., 
that the source was detected regardless of the number of 
principal components). We then listed sources with a signal-to- 
noise ratio (S/N) greater than 5 (or with a signal exceeding 4c) 


and determined the optimal number of principal components 
(comp) that maximized the S /N for each source. 


3.3. Source Magnitude Calculation 


To calculate the apparent and absolute magnitudes of the 
detected sources, we proceeded as follows: (1) We fit a 2D 
Gaussian profile to the point-spread function (PSF) cube to 
obtain the FWHM in milliarcseconds (mas) and to recenter the 
PSF frames. (2) The PSF cube was reduced into a single 2D 
PSF by computing the median of the frames. (3) We 
normalized the PSF so that the flux within a radius of 1 
FWHM equated to 1. (4) Given a thermal artifact affecting the 
quality of the real PSF in our 2020 observations (as discussed 
in Prasow-Émond et al. 2022), we generated a normalized 
synthetic 2D PSF using the FWHM. (5) From the list of 
sources with S/N > 5, we generated a list of approximate 
coordinates using ds9. (6) Using the preprocessed cube, the 
optimal number of principal components, and the approximate 
coordinates, we fit the astrometric parameters (0, the position 
angle in degrees, and p, the relative separation from the X-ray 
binary in pixels) and photometric parameters (fi, the number of 
counts within an aperture radius of 3 FWHM, i.e., the relative 
flux) using a Nelder—Mead optimization algorithm from VIP. 
This fitting process, involving the injection of synthetic sources 
with a negative flux at source location (Lagrange et al. 2010), 
aimed to minimize X? residuals. (7) Once the parameters were 
determined, the position angle in the images was converted into 
the true position angle using the celestial north of NIRC2, 
which is 02262 + 02018 (Service et al. 2016; see Franson et al. 
2022 for an example and more details). (8) We converted the 
units of p from pixels to mas using the plate scale (see 
Section 3.1). (9) To derive the apparent magnitude of the 
sources (Mec), we applied the following equation: 


Mec = —2.5 logig( f/f) + mxns. (1) 


where f» is the PSF flux within the same aperture as the cube (3 
FWHM) and mygg is the apparent magnitude of the X-ray 
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Table 1 
Information on the Binary Nature for the 14 Observed X-Ray Binaries 
Target Type Class Donor Star References Compact Object References 
(1) (2) (3) (5) (6) (7) 
RX J1744.7—2713 HMXB N/A BO.5III-Ve (1) WD* (2) 
BeXRB 
y Cas analog 
IGR J18483—0311 HMXB Outburst B0.5-BI/BO-1Iab (3, 4) NS (5) 
SEXT 
XP 
y Cas HMXB Steady BO.5IVe (6, 7) WD* (8) 
BeXRB 
SAX J1818.6—1703 HMXB Variable BO.5Iab (4) NS* (9) 
SFXT 
1H 2202+501 HMXB N/A B3Ve (10) N/A 
BeXRB* 
4U 2206+543 HMXB Variable O9.5V (11) NS (12) 
4U 1700—377 sgHMXB Variable O6lafpe (13) NS (14) 
IGR J17544—2619 HMXB Flaring O9Ib/O9IV-V (15, 16, 17) NS (18) 
SEXT 
RX J2030.5+4751 HMXB N/A BO.5III-Ve (19) N/A 
BeXRB 
y Cas analog 
Cyg X-1 HMXB Variable O9.7Iabpvar (20) BH Q1) 
Microquasar 
X Per HMXB Variable BOVe (22) NS (23) 
BeXRB 
XP 
1H 0556+286 HMXB* N/A B5ne (24) WD or NS (25) 
RX J0648.1—4419 HMXB N/A sdO5.5 (26) WD (27) 
XP 
Vela X-1 HMXB Variable BO.5Ib (28) NS (28) 
XP 


Notes. Column (1): name of the target. Column (2): subclass /type of the X-ray binary, as found in Bird et al. (2016). Column (3): class of X-ray emission, as found in 
Krimm et al. (2013). Column (4): the spectral type of the donor star. Column (5): reference for the donor star. Column (6): the nature of the compact object. Column 


(7): reference for the compact object. 
* Indicates an uncertain nature. 


References. (1) Sarty et al. 2011; (2) Lopes de Oliveira et al. 2006; (3) Chaty et al. 2008; (4) Torrejón et al. 2010; (5) Sguera et al. 2007; (6) Moffat et al. 1973; (7) 
Raguzova & Popov 2005; (8) Postnov et al. 2017; (9) Walter & Zurita Heras 2007; (10) Simon et al. 2019; (11) Negueruela & Reig 2001; (12) Torrejón et al. 2004; 
(13) Sota et al. 2014; (14) Reynolds et al. 1999; (15) Pellizza et al. 2006; (16) Giménez-Garcia et al. 2016; (17) Bikmaev et al. 2017; (18) in’t Zand 2005; (19) Motch 
et al. 1997; (20) Sota et al. 2011; (21) Bolton 1972a; (22) Lyubimkov et al. 1997; (23) White et al. 1977; (24) Popper 1950; (25) Liu et al. 2006; (26) Jaschek & 


Jaschek 1963; (27) Popov et al. 2018; (28) Hiltner et al. 1972. 


binary (see Section 3.3.2). (10) Using the known distance in 
parsecs (see Section 3.3.1), the apparent magnitude of the 
candidate CBCs (me) was converted to absolute magnitude 
(Mecc) via the distance modulus equation. (11) Finally, to 
estimate the mass of the sources, we compared the absolute 
magnitude (Mec) with evolutionary models from MESA 
Isochrones & Stellar Tracks (MIST; Paxton et al. 2011, 2013, 
2015, 2018; Choi et al. 2016; Dotter 2016; Faherty et al. 2016) 
at the system’s age (see Section 3.3.3). 

Due to the poor quality of the PSF in our 2020 observations, 
we used the synthetic 2D PSF for fitting processes related to 
these observations. However, we kept the original PSF for our 
2017 and 2018 observations. We conducted tests using both 


synthetic and real PSFs on the 2017 and 2018 data, yielding 
consistent results. Consequently, the use of a synthetic PSF 
does not affect significantly the derived parameters. 

The upcoming sections detail the acquisition of parameters 
used in the magnitude calculations. 


3.3.1. Determining the Distance from the Observer 


The distance from the observer is presented in the second 
column of Table 5 in the Appendix. Distances for 1H 2202 
+501, 4U 22064-543, 4U 1700—377, IGR J17544—2619, Cyg 
X-1, X Per, and Vela X-1 were obtained from Zhao et al. 
(2023), which uses parallax measurements to infer distances 
using either an inversion or a Bayesian approach for a catalog 
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Table 2 
Keck/NIRC2 Observing Log 
UT Date Target Filter WFS tint Coadds Neames P.A. Cov. 
(s) (deg) 
2017 Sep 8 RX J1744.7—2713 L’ SH 0.5 60 40 14.3 
Cygnus X-1 L' SH 0.5 60 140 60.6 
y Cassiopeiae L’ SH 0.18 150 130 61.2 
X Persei L' SH 0.5 60 89 43.9 
2018 Jan 3 X Persei L! SH 1 45 54 36.7 
1H 0556+286 L’ SH 1 45 50 48.6 
RX J0648.1—4419 L' SH 1 45 17 6.5 
Vela X-1 Et SH 1 45 60 26.0 
2020 Jul 11 RX J1744.7—2713 L' py 0.5 60 120 38.3 
IGR J18483—0311 L' py 0.4 60 125 42.6 
y Cassiopeiae L' py 0.0528 400 150 48.1 
2020 Jul 12 RX J1744.7—2713 K, py 0.6 45 92 39.6 
SAX J1818.6—1703 L' py 0.4 60 90 26.3 
1H 2202+501 BS py 0.4 60 89 36.9 
4U 2206+543 L py 0.4 60 27 49.5 
2020 Jul 13 4U 1700—37 L SH 0.4 50 69 15.1 
IGR J17544—2619 L' SH 0.4 50 94 36.5 
RX J2030.5+4751 E SH 0.4 60 132 53.9 
4U 2206+543 L' SH 0.4 60 149 49.5 


Note. Abbreviations are as follows: WFS—wave front sensor; SH—Shack-Hartmann; py—pyramid WFS; 1;,,—integration time; Ngames—number of frames; P.A. 


Cov.—parallactic angle coverage. 


of X-ray binaries. For the other X-ray binaries, the distance was 
estimated through a photogeometric calculation (Bailer-Jones 
et al. 2021). It was calculated using the parallax measurement 
and its uncertainty (geometric), as well as the G magnitude and 
the BP — RP color (photometric) from the Gaia Data Release 3 
(DR3; De Angeli et al. 2023; Montegriffo et al. 2023). Table 3 
presents the Gaia DR3 ID for each target. In cases where an 
object's distance was sourced from the literature, the respective 
reference is cited in the literature review of Section 2. As 
previously mentioned, this study targets X-ray binaries within 
~2-3 kpc accessible with Keck/NIRC2. 


3.3.2. Determining the Apparent Magnitude of the Central X-Ray 
Binary 


Observing X-ray binaries in the L’ band of Keck/NIRC2 is 
not standard practice, making the direct determination of the 
true apparent magnitudes of the central X-ray binaries ("mxgg) 
unfeasible. Nonetheless, we leveraged the W1 filter from the 
Wide-field Infrared Survey Explorer (WISE), which has a 
central wavelength similar to Keck/NIRC2 L’ band 
(A = 3.353 um). Using the WISE Source Catalog (Wright 
et al. 2010), we approximated mxgg for all X-ray binaries. 
These values are presented in the fifth column of Table 5 in the 
Appendix. To account for the differences between the two 
filters, we considered an uncertainty of 0.5 mag on rxgg. 


3.3.3. Determining the Age of the System 


Most X-ray binaries in our sample lack age estimates in the 
literature, except for RX J1744.7—2713 (up to ^60 Myr; 
Coleiro & Chaty 2013), 4U 1700—377 (up to ^80 Myr; 
Coleiro & Chaty 2013), y Cas (8.0 + 0.4 Myr; Zorec et al. 
2005), Cyg X-1 («4 Myr; Miller-Jones et al. 2021), and X Per 
(~5 Myr; Lyubimkov et al. 1997). For the remaining X-ray 


Table 3 
Gaia DR3 ID for Each Target 

Target Gaia DR3 ID 

RX J1744.7—2713 4060784345959549184 
IGR J18483—0311 4258428501693172736 
y Cas 426558460884582016 
SAX J1818.6—1703 4097365235226829312 
1H 22024-501 1979911002134040960 
4U 2206--543 2005653524280214400 
4U 1700—377 59763829 15813535232 
IGR J17544—2619 40639088 10076415872 
RX J2030.5+4751 2083644392294059520 
Cyg X-1 2059383668236814720 
X Per 168450545792009600 
1H 05564-286 3431561565357225088 
RX J0648.1—4419 5562023884304070000 
Vela X-1 5620657678322625920 


binaries, we established an upper limit using basic evolutionary 
models and the spectral type of the donor star (see Table 1). We 
found the maximum age that the donor star can reach before 
exploding in supernovae. However, this approach provides 
only an approximate estimation of the age. These values are 
presented in the sixth column of Table 5 in the Appendix. 


3.3.4. Determining the Errors 


Errors in the fit parameters (0, p, and fı) were estimated 
using an injection/recovery approach (see Prasow-Émond et al. 
2022). This approach relies on injecting synthetic sources with 
known parameters into the images. Subsequently, the same 
optimization method was applied (see Section 3.3), and the 
error was determined as the difference between estimated and 
known parameters. This method was employed across a range 
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of parameter values to ensure consistency; it was observed that 
errors were more pronounced for smaller offset values (i.e., 
those closer to the central X-ray binary). Additional sources of 
uncertainty were taken into consideration in cases where the 
value of a parameter is expected to remain consistent between 
two sets of observations, specifically astrometric parameters in 
two different bands. The dominant source of uncertainty was 
defined as the total uncertainty. 


4. Key Properties of the Sample 
4.1. ^; Cassiopeiae 


y Cassiopeiae—also known as 2S 0053+604 (hereafter ^ 
Cas)—harbors a well-studied optical component classified as a 
Be star (Moffat et al. 1973). Its X-ray luminosity (—10??— 
l0? ergs !; Raguzova & Popov 2005) is higher than the 
typical luminosity for O/B stars (~10°° ergs ') but too low to 
be a Be/NS system (~10*4 erg s^! in quiescence; Shrader et al. 
2015). The nature of the system can be explained by two 
hypotheses: (1) the system is an HMXB, involving accretion 
onto a WD or a fast-spinning NS (Postnov et al. 2017); or (2) 
the excess of X-ray emission stems from physical processes in 
the high atmosphere of y Cas (Kubo et al. 1998; Robinson & 
Smith 2000). Though y Cas is not confirmed as being an 
HMXB, its resemblance to other sources, referred to as y Cas 
analogs, warranted its inclusion in our sample. Located at a 
distance of 0.19+0.02 kpc, it is the nearest system in our 
sample. Its proximity and proper motions (25.7 + 0.5 mas yr! 
in R.A., —3.9+0.4 mas yr ' in decl.; Perryman et al. 1997) 
allowed us to conduct a proper-motion analysis within the 
interval between our two observation sets (see Section 6.2). 


4.2. RX J1744.7-2713 


RX J1744.7-2713 is classified as a BeXRB (Israel et al. 
1997) and is composed of a BO.5III- Ve star (Motch et al. 1997; 
Steele et al. 1999; Lopes de Oliveira et al. 2006) and a WD. 
However, the origin of X-ray emission is still uncertain and 
debated in the literature (Lopes de Oliveira et al. 2006). It is 
also known as a y Cas analog, due to the similarities in their 
X-ray properties (e.g., Shrader et al. 2015). The first high- 
contrast images of this HMXB were presented in Prasow- 
Émond et al. (2022), in which more comprehensive information 
on the system can be found. 


4.3. 4U 1700—377 


4U 1700—377, discovered with the Uhuru X-ray satellite, is 
classified as an HMXB (Jones et al. 1973). The system contains 
a supergiant donor star of spectral type O6lafpe (Sota et al. 
2014) and a magnetized NS (e.g., Reynolds et al. 1999; Bala 
et al. 2020; van der Meij et al. 2021) exhibiting strong flaring 
activity (e.g., Kuulkers et al. 2007). This HMXB was observed 
with, e.g., XMM-Newton (van der Meer et al. 2005; Giménez- 
García et al. 2015), Chandra (Boroson et al. 2003; Martínez- 
Chicharro et al. 2021), the Hubble Space Telescope (Hainich 
et al 2020) the Fiber-fed Extended Range Optical 
Spectrograph (Hainich et al. 2020), IUE (Dupree et al. 1978), 
EXOSAT (Haberl et al. 1989), BATSE (Rubin et al. 1996), and 
BeppoSAX (Reynolds et al. 1999). NGC 6231, located within 
the OB association Sco OBI, has been recently confirmed as 
the parent cluster of 4U 1700—377 (van der Meij et al. 2021). 
Moreover, Coleiro & Chaty (2013) inferred an estimated age of 
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~80 Myr for the system, a value we have regarded as an upper 
limit as in Prasow-Emond et al. (2022). Its X-ray luminosity 
can reach up to ~7 x 10°% erg s^! (Laurent et al. 1992). 


4.4. 4U 2206+543 


Although the nature of its donor star remains uncertain (first 
believed to be of spectral type Be, Steiner et al. 1984; then 
O9.5V, Negueruela & Reig 2001), 4U 2206+543 is a well- 
studied HMXB harboring a magnetar NS (e.g., Torrejón et al. 
2004, 2018; Blay et al. 2005; Reig et al. 2009; Ikhsanov & 
Beskrovnaya 2010). First discovered with Uhuru (Giacconi 
et al. 1972), the system was subsequently observed with, e.g., 
EXOSAT (Saraswat & Apparao 1992), RXTE (e.g., Corbet & 
Peele 2001), IBIS/ISGRI on INTEGRAL (e.g., Bird et al. 
2004), BeppoSAX (e.g., Masetti et al. 2004), the Very Large 
Array (Blay et al. 2005), Swift (Corbet et al. 2007), and Suzaku 
(Finger et al. 2010). Its X-ray luminosity ranges from 
~10°% erg s^! in quiescence up to ~10*°-10*° erg s! during 
more active phases (Ribó et al. 2006). 


4.5. RX J2030.5--4751 


RX J2030.5+4751, also known as BD +47°3129 and SAO 
49725, was first discovered with ROSAT (Motch et al. 1997). It 
is identified as a BeXRB and y Cas analog. It has a BO.SIII-Ve 
spectral type donor star and a maximum X-ray luminosity 
of ~10%° ergs ' (Liu et al. 2006; Raguzova 2007). XMM- 
Newton/EPIC observations have revealed that RX J2030.5 
+4751 has a hard X-ray spectrum, suggesting the presence of a 
dense, large, and stable circumstellar environment surrounding 
the system (Lopes de Oliveira et al. 2006). On 2016 July 20, RX 
J2030.5+4751 underwent type I (i.e., smaller and repetitive) 
bursts, reaching its maximum luminosity (progressive weaken- 
ing since; Steele 2016a, 2016b, 2016c). Regarding the long-term 
variability, Servillat et al. (2012) reported a significant 
nonperiodic variability of approximately 1 mag in the light 
curve of RX J2030.5+4751 over ~100 yr, likely caused by 
changes in the properties of the decretion disk. 


4.6. 1H 2202+501 


1H 22024-501 is a poorly studied HMXB, albeit appearing 
in some surveys (e.g., Wood et al. 1984; Hanson et al. 1996; 
Liu et al. 2000). It consists of a Be star of spectral type B3Ve 
(Simon et al. 2019), and the nature of the compact object 
remains uncertain. Its X-ray luminosity is estimated to be 
^9 x l0? erg s! (Chevalier & Ilovaisky 1998). 


4.7. SAX J1818.6—1703 


SAX J1818.6—1703 was discovered with BeppoSAX while 
undergoing a strong outburst (in't Zand et al. 1998). 
Subsequent observations were carried out using IBIS/ISGRI 
on board INTEGRAL (e.g., Grebenev & Sunyaev 2005; Sguera 
et al. 2005; Zurita Heras & Chaty 2009; Sidoli et al. 2016), 
RossiXTE (Sguera et al. 2005; Smith et al. 2012), and Swift 
(e.g., Bird et al. 2009; Sidoli et al. 2009). The system has 
similar properties to other SFXTs (e.g., Sguera et al. 2005; 
Negueruela & Smith 2006; Sidoli et al. 2009; Bozzo et al. 
2012), and its compact object is likely an NS (e.g., Walter & 
Zurita Heras 2007; Boon et al. 2016). The donor star of the 
system, confirmed by Chandra (in't Zand et al. 2006), is 
classified as a supergiant star of spectral type BO.5Iab (Torrejón 
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et al. 2010). Also observed using XMM-Newton, SAX J1818.6 
—1703 has a quiescent X-ray luminosity that can drop to values 
below ~10* erg s^! (as determined by not being detected in 
Bozzo et al. 2012) and can reach up to ~8 x 10%” ergs | 
(Torrejón et al. 2010). 


4.6. IGR J18483 —0311 


Discovered with IBIS/ISGRI on board INTEGRAL by 
Chernyakova et al. (2003), IGR J18483—0311 is a well-studied 
SEXT (Rahoui & Chaty 2008) composed of an X-ray pulsar 
(Sguera et al. 2007) and a supergiant donor star of spectral type 
BO.5Ia (Chaty et al. 2008) or B0.5-B 1Iab (Torrejón et al. 
2010). The system undergoes multiple short and long outbursts 
(e.g., Sguera et al. 2007, 2010; Ducci et al. 2013; Sguera et al. 
2015), resulting in its X-ray luminosity ranging from ~10°*— 
10°* erg s^! in quiescence (Romano et al. 2010; Sguera et al. 
2015) up to ~10*° erg s^! during its strongest flares (Torrejón 
et al. 2010). By considering evolutionary scenarios and 
exploring the relationship between spin and orbital periods, 
Liu et al. (2011) suggested that the compact object in IGR 
J18483—0311 may originate from an O-type emission-line star 
rather than an average main-sequence star. 


4.9. IGR J17544 —2619 


IGR J17544—2619 was first discovered near the Galactic 
center while undergoing short (a few hours) outbursts using 
IBIS/ISGRI on board INTEGRAL (Grebenev et al. 2003, 
2004; Sunyaev et al. 2003). The system was subsequently 
observed with, e.g., XMM-Newton (Drave et al. 2014; 
González-Riestra et al. 2004), Chandra (in't Zand 2005), 
EMMI/SOFI/NTT (Pellizza et al. 2006), Suzaku (Rampy et al. 
2009), RXTE (Drave et al. 2012), Swift (e.g., Romano et al. 
2015), and NuSTAR (e.g., Bhalerao et al. 2015). It is identified 
as an SFXT (Negueruela et al. 2006), and the optical/near- 
infrared counterpart is a massive star of spectral type O9Ib 
(25-28 Mo; Pellizza et al. 2006; Giménez-García et al. 2016) or 
O9IV-V (23 Mo; Bikmaev et al. 2017). The compact object is 
an NS (in't Zand 2005), as inferred by the presence of a 
cyclotron line at 17keV and the magnetic field strength 
(~1.5 x 10? G, typical for NSs in X-ray binaries; Bhalerao 
et al. 2015). The change in X-ray luminosity between 
quiescence and outburst is significant, ranging from 
Lx ~ 10°°-10** erg s! (e.g., Bozzo et al. 2016) during 
quiescence to a maximum of Ly-3 x 10°% ergs ! during 
outburst (Romano et al. 2015). 


4.10. Cyg X-1 


Cygnus X-1 (hereafter Cyg X-1), first discovered in 1964, is 
one of the most well-studied astronomical objects (e.g., Fabian 
et al. 1989; Esin et al. 1998; Orosz et al. 2011; Tomsick et al. 
2014; Parker et al. 2015; Sell et al. 2015; Walton et al. 2016; 
Mastroserio et al. 2019). The system's compact object, with a 
mass of 21.2 + 2.2 Mọ (Miller-Jones et al. 2021), was the first 
observed candidate BH (Murdin & Webster 1971; Bolton 
19722), leading to significant breakthroughs in the astronom- 
ical scientific community. Cyg X-1 is classified as an HMXB 
(Bolton 1972b), and the donor star is characterized as a massive 
star of spectral type O9.7Iabpvar (Sota et al. 2011). The system 
is fairly young, with an estimated age of 5+ 1.5 Myr in 
Mirabel & Rodrigues (2003) and «4 Myr in Miller-Jones et al. 
(2021). Cyg X-1 is located at FIIN kpc; its distance from 
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the observer was precisely estimated in Miller-Jones et al. 
(2021) based on radio parallax measurements and validated 
with Gaia DR3. The microquasar undergoes variable X-ray 
emission (Krimm et al. 2013), with a maximum X-ray 
luminosity of ~3 x 10°” erg s^! (Di Salvo et al. 2001). 


4.11. Vela X-1 


Vela X-1 is a pulsing HMXB discovered with the Uhuru 
X-ray satellite (Giacconi et al. 1971) and was observed through 
multiple surveys and with several instruments (e.g., Charles 
et al. 1978; La Barbera et al. 2003; Fürst et al. 2014; Martínez- 
Nüfiez et al. 2014). The system is highly variable and 
undergoes transient outbursts and X-ray eclipses (e.g., van 
der Klis & Bonnet-Bidaud 1984). Its X-ray luminosity can be 
as high as 4 x 10°% ergs ^! (Kreykenbohm et al. 2008). It is 
composed of a B0.51b donor star (Hiltner et al. 1972) and of a 
magnetized NS (e.g., Hiltner et al. 1972; van Kerkwijk et al. 
1995; Kreykenbohm et al. 2002; Diez et al. 2022). A complete 
review of this object can be found in Kretschmar et al. (2021). 


4.12. RX J0648.1 —4419 


RX J0648.1—4419 is a unique X-ray pulsating system, as it 
is the only HMXB known to date containing a hot subdwarf 
donor star of spectral type O (sdO5.5; Jaschek & Jaschek 1963; 
Mereghetti et al. 2009). The compact object has a mass of 
1.28 +0.05 Mo (Mereghetti et al. 2009) and was initially 
believed to be an NS (e.g., Israel et al. 1997; Mereghetti et al. 
2016). Popov et al. (2018) suggested that it was likely a young 
(~2 Myr) contracting WD. Over the course of almost 30 yr of 
observations (e.g., Mereghetti et al. 2011, 2021; La Palombara 
et al. 2015), the X-ray luminosity remained stable, maintaining 
a value of ^10? ergs |. 


4.13. 1H 0556+286 


1H 0556+286 contains a Be star of spectral type B5ne 
(Popper 1950). It is a poorly studied system; while it is 
generally thought to be an HMXB (e.g., Helfand & Moran 
2001; Liu et al. 2006), Torrejón & Orr (2001) presented 
BeppoSAX observations in which no X-ray emission was 
detected. As per these results, neither the Be/NS nor Be/WD 
scenarios appear probable; hence, the nature of the system 
remains unclear. 


4.14. X Per 


X Persei (hereafter X Per), also known as 4U 0352+309, 
was discovered in 1972 with the Copernicus Observatory 
(Hawkins et al. 1975; Mason et al. 1976). The system was 
subsequently observed with, e.g., the High Energy Astro- 
physical Observatory (HEAO 2/Einstein; Weisskopf et al. 
1984), RXTE (e.g., Delgado-Martí et al. 2001; Coburn et al. 
2002), INTEGRAL (e.g., Lutovinov et al. 2012), XMM- 
Newton (e.g., La Palombara & Mereghetti 2007), and Chandra 
(e.g., Valencic & Smith 2013). The system is identified as an 
HMXB/BeXRB, composed of a magnetized NS as the 
compact object (e.g., White et al. 1977; Coburn et al. 2001; 
Doroshenko et al. 2012; Maitra et al. 2017; Yatabe et al. 2018) 
and a Be star of spectral type BOVe as the donor star 
(Lyubimkov et al. 1997). In quiescence, its X-ray luminosity is 
~10* erg s^! (e.g., Coburn et al. 2001). When undergoing 
strong outburst activity, its X-ray luminosity can reach up to 
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~2 x 10°% erg s ! (Lutovinov et al. 2012). Similarly to Cyg 
X-1, the system is also relatively young, with an estimated age 
of ~5 Myr (Lyubimkov et al. 1997). Objects with similar 
properties are often referred to as X Per analogs. 


5. High-contrast Images 
5.1. ^j Cas 


As discussed in Section 4.1, observations of y Cas in the L’ 
band were initially made on 2017 September 8 and 
subsequently revisited almost 3 yr later on 2020 July 11. 
Figure 2 presents the L’-band high-contrast images of ~y Cas for 
both epochs. A bright source, labeled B, was detected with an 
S/N > 5. A much fainter source, labeled C, was also detected, 
with an S/N ~ 3. In Section 6.2, we undertake a proper-motion 
analysis to determine whether these sources are more likely to 
be bound CBCs or background stars. 


5.2. Other X-Ray Binaries 


Figure 3 presents a panel of L'-band high-contrast images of 
all the other X-ray binaries, in order: 4U 1700—377, 4U 2206 
+543, RX J2030.5+4751, 1H 2202+501, SAX J1818.6 
—1703, IGR J18483—0311, IGR J17544—2619, Cyg X-l, 
Vela X-1, RX J0648.1—4419, 1H 0556+286, and X Per. The 
L'- and K,-band images of RX J1744.7—2713 can be found in 
Prasow-Emond et al. (2022). By inspecting the images, we 
determined that 4U 2206+543, Vela X-1, RX J0648.1—4419, 
and 1H 0556+286 do not exhibit any potential candidate 
CBCs. Consequently, further analysis of these systems will not 
be pursued, except for the calculation of companion frequency 
in Section 6.3. Among the remaining X-ray binaries, we 
successfully detected at least one source for each system with a 
significantly large S/N (>5). These sources are labeled in the 
images, starting from the letter B. 

Table 6 in the Appendix presents several physical properties 
of the detected sources, including the angular separation in 
mas, the position angle in degrees, the apparent magnitude in 
L’ band, and the mass estimated from evolutionary models. 

In the next sections, we analyze the nature of the detected 
sources and discuss the implications of the results. 


6. On the Nature of the Detected Sources 
6.1. Background Contamination 


We used TRILEGAL (Girardi et al. 2005) to discriminate 
background/foreground stars from gravitationally bound 
CBCs. TRILEGAL is a 3D model employed to simulate the 
photometric properties of star fields within the Galaxy (e.g., 
Chun et al. 2015; Dietrich & Ginski 2018; Jones et al. 2021; 
Williams et al. 2021). We compiled a list of predicted sources 
within a 1 x 1 arcmin? region surrounding the X-ray binary 
using its R.A./decl. coordinates. Subsequently, we applied a 
geometric rescaling so that the TRILEGAL FOV matches the 
FOV of the high-contrast images (10.24 x 10.24 arcsec’) for 
Statistical consistency. We calculated the cumulative distribu- 
tion of the expected number of sources in the FOV below an 
apparent magnitude (mz), for the distribution of apparent 
magnitudes (L’). It is denoted as ngoy(L’ < my’). 

Figure 4 presents ngoy(L’ < my’) as a function of mz for 
both TRILEGAL and the detected sources in the high-contrast 
images. It includes all X-ray binaries with detected sources, 
except Cas (for which the confirmation of the nature of the 
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CBCs relies primarily on proper-motion analysis) and RX 
J1744.7—2713 (presented in Prasow-Emond et al. 2022). In 
Appendix Table 6, we listed npoy(L! € my) source) for each 
detected source. This represents the expected number of 
sources—based on TRILEGAL simulations—with apparent 
magnitudes (mz) below the magnitude of the corresponding 
source (ML) source) While accounting for errors. If the calculated 
value, including the upper limit, was lower than the total 
number of sources detected below mz) source, We reject the 
hypothesis of background/foreground contamination for that 
particular source. 

Many sources were not predicted by TRILEGAL. Specifi- 
cally, all detected sources in 4U 1700—377, RX J2030.5 
+4751, Cygnus X-1, X Per, and 1H 2202+501 were not 
expected from the model given their magnitudes and would 
thus suggest that they are bound to the X-ray binary. The FOVs 
of IGR J17544—2619 and SAX J1818.6—1703 were expected 
to be more populated than what we detected, suggesting that 
the detected sources are likely background or foreground 
contaminants. The disparity between the predicted and 
observed number of sources might result from the elimination 
of stationary sources during the ADI process or from an 
insufficient S/N for detection. In the case of IGR 718483 
—0311, sources B, C, and F might be CBC candidates, but their 
status remains uncertain in the current stage of our study. 

We extended our analysis by calculating the probability of 
chance alignment for each detected source. This probability 
represents the likelihood that these sources are not associated 
with the X-ray binary system based on the angular separation 
and the density on the sky of unrelated objects. This method 
assumes that the distribution of unrelated sources across the 
area follows a Poisson distribution. Note that this method 
usually relies on sky surveys such as the Two Micron All Sky 
Survey Point Source Catalog (e.g. Correia et al. 2006; 
Lafreniére et al. 2008, 2014; Prasow-Emond et al. 2022) to 
establish the distribution of unrelated objects. However, in this 
work we used TRILEGAL as an alternative owing to the 
absence of available K,-band observations. 

To calculate the aforementioned probability, we first 
divided the cumulative distribution of the number of sources 
NtTRILEGAL(L’ < mr) by the area from which the sources were 
retrieved (between 6x6 arcmin? and 15 x 15 arcmin? 
depending on the location of the X-ray binary). This division 
enabled us to derive a surface density denoted as X. Using the 
angular separation O in arcsec, the probability of a source 
being drawn from the TRILEGAL distribution—thus indicat- 
ing its lack of association with the central X-ray binary—is 
given by 


Pynrelated 27, ©) =1- exp( —5530?). (2) 


In Appendix Table 6, we listed 1 — Pynrelatea(, ©) as 
percentages. Many sources have high probabilities (78596) of 
being associated with the central X-ray binary. However, some 
sources have lower probabilities (between 65% and 75%) that 
are not as statistically significant, but we nevertheless identified 
them as candidate CBCs given the early stage of the study. 
Sources with probabilities below 60% (not statistically 
significant) were rejected as candidate CBCs. 


6.2. Proper-motion Analysis: \ Cas 


Conducting a proper-motion analysis is among the most 
robust methods for confirming the gravitational association 
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Keck/NIRC2 L'-band 
2017 Sept 08 
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Keck/NIRC2 L'-band 
2020 July 11 


0.5" = 95 AU 


Figure 2. Keck/NIRC2 L'-band high-contrast images of y Cas acquired on 
2017 September 8 (up) and 2020 July 11 (bottom), treated and reduced using a 
PCA annular ADI algorithm (using VIP; Gomez Gonzalez et al. 2017). The 
sources detected with S/N > 3 are labeled B and C. The white cross denotes 
the approximate position of y Cas, masked by the coronagraph. The insets in 
the lower left corners show zoomed-in high-contrast images with a focus on an 
annular region (obtained using the PCA annulus algorithm in VIP), 
highlighting the presence of the fainter source (labeled as C). 


of a source with the system (e.g., Bohn et al. 2021). Given the 
proximity of y Cas and its proper motions (see Section 4.1), 
it was possible to conduct a statistically significant proper- 
motion analysis between the two observed epochs (2017 
September 8 and 2020 July 11). Figure 5 presents the relative 
separations between sources B, C, and y Cas in R.A. and 
decl, alongside the expected position of a stationary 
background star. The figure also displays the angular 
separation and position angle over time, along with the 
expected tracks for both comoving and stationary background 
objects. 
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Source B's trajectory implies an underlying motion that 
necessitates additional epochs of observation for validation. In 
2020, the angular separation data point aligns with the comoving 
track, but the position angle data point deviates by ~3o from the 
same track. This trajectory suggests that source B is more likely 
to be bound to y Cas rather than an unrelated background or 
foreground object. Nonetheless, its motion suggests potential 
scenarios such as orbital motion, ejection from the system, or the 
presence of a nonstationary background or foreground object. 
Using the mass of y Cas (~13 Mo; Nemravová et al. 2012) and 
the radial separation of source B, we calculated the escape 
velocity as v. = /2GM/p ~ 7680 ms". Additionally, the 
projected velocity between 2017 and 2020 was determined, 
resulting in the vector Vproj zz (1670 p — 17964 0) m s |. Since 
the norm of the velocity vector is greater than the escape 
velocity, source B appears to be physically associated with ^ 
Cas, but it is not bound (as also suggested in Hutter et al. 2021). 
A comprehensive characterization of this motion using high- 
contrast imaging necessitates further epochs of observation. 

Source C consistently follows the motion track of an 
unrelated object across all three plots. Thus, we excluded it 
from our list of candidate CBCs with a >>3o confidence level. 


6.3. Frequency of Multiple Systems and Companion Frequency 


Two key concepts are commonly defined in the stellar 
multiplicity literature (e.g., Duchéne & Kraus 2013): the 
frequency of multiple systems (MF) and the companion 
frequency (CF; i.e., the average number of companions per 
target). While a proper-motion analysis is required to confirm 
most of the sources, we calculated a first estimation of MF and 
CF for the observed X-ray binaries in our sample. Among the 
total of 14 observed X-ray binaries, we have identified 
candidate CBCs in eight systems: 4U 1700—377, RX 
J2030.54-4751, Cyg X-1, X Per, 1H 2202+501, y Cas, RX 
J1744.7—2713, and IGR J18483—0311. Thus, based on these 
numbers and at this stage of the study, MF for triple or higher- 
tier systems would be 8/14~0.6 (~60%). For CF, it is 
important to note that X-ray binaries inherently possess one 
companion, the donor star, which is included in our calculation 
of CF. Based on the status of the detected sources listed in 
Appendix Table 6, the calculated average number of 
companions per compact object is 2.1 + 1.1 (210% + 110%). 
This means that, on average, the compact object has two 
companions (the donor star and one additional companion). 
However, this value reduces to 1.8 + 0.9 (180% + 90%) when 
considering only the sources that are most likely to be 
gravitationally bound (1 — Punrelatea( ©, ©) 8596). These 
values are subject to change as new observations become 
available and further analyses are conducted. 


7. Discussion 
7.1. Stellar Multiplicity 


The discovery of candidate CBCs would imply that X-ray 
binaries can still be produced by multiple-star systems rather 
than exclusively binary systems. The total mass (i.e., compact 
object and donor star) of all the HMXBs in our sample exceeds 
^10 Mọ, with some reaching up to around 60 Mọ, placing 
them on the higher end of the mass spectrum. Stellar 
multiplicity is believed to be common in high-mass star 
systems (e.g., Chini et al. 2012; Duchéne & Kraus 2013), and 
high-order multiplicity is thought to increase with primary 
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Figure 3. Keck/NIRC2 L’-band high-contrast images of all observed X-ray binaries—except RX J1744.7—2713 (see Prasow-Emond et al. ) and y Cas (see 
Figure 2)—acquired on 2017 September 8, 2018 January 3, and 2020 July 11-13. Images were treated and reduced using a PCA annular ADI algorithm (using VIP; 
Gomez Gonzalez et al. 2017). The sources detected with S/N > 5 when computing the S/N map are labeled. The white cross denotes the position of the X-ray binary 
masked by the coronagraph. North points upward, and east points to the left, as in Figure 
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Figure 4. Cumulative distribution of the number of sources in the FOV expected from TRILEGAL simulations (red) and detected with an S/N > 5 in the high- 
contrast images (green) as a function of the apparent magnitude mz. It includes all X-ray binaries in which sources were detected except y Cas and RX J1744.7—2713 
(see Prasow-Emond et al. 2022): 4U 1700—377, RX J2030.5+4751, IGR J17544-2619, Cygnus X-1, 1H 2202+501, IGR J18483—0311, X Per, and SAX 


J1818.6—1703. 


mass (e.g., Peter et al. 2012). However, surveys for high-mass 


stars remain incomplete. 


For high-mass stars (716 Mc), MF and CF are estimated to 
be 28096 and 130% + 20%, respectively (Chini et al. 2012; 
Duchéne & Kraus 2013). The first estimation of MF for our 


sample (~60%; calculated in Section 6.3) falls below this 
percentage. However, in this study, MF is constrained by the 
range of projected separations (up to ~12,000 au). This implies 
that increasing this limit could potentially lead to the discovery 
of more companions and hence increasing the estimation of 
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Figure 5. Left: relative separations between source B (top row), source C (bottom row), and y Cas in R.A. (o) and decl. (6). The first-epoch astrometric point is plotted 
in blue (2017 September 8), and the second-epoch astrometric point is plotted in red (2020 July 11). The expected position for a stationary background object is plotted 
in yellow, along with its proper-motion track. Middle: separation from y Cas in mas as a function of time. A background object with zero proper motion would follow 
the orange track, while a bounded and comoving CBC would lie within the blue zone. Right: same as the middle panels, but for the position angle in degrees. 


MF. As for CF (210% + 110%; calculated in Section 6.3), it 
currently exceeds the estimate obtained from the literature 
(130% + 20%; Chini et al. 2012). Further observations will 
likely lead to the rejection of candidate CBCs we have 
detected, which would lower the sample's CF (along with the 
associated uncertainty range). In both cases, our preliminary 
estimations (MF and CF) seem to be broadly in line with 
current estimates in the literature. 

For solar-type stars, the frequency N(n) of multiplicity n 
follows a geometric distribution N(n) ~ 8 " (up to n— 7 with 
B — 2.3 or 3.4; Eggleton & Tokovinin 2008). Such a relation has 
not been derived for massive stars. Thus, the results from this 
study could be used to derive one, allowing us to predict the 
frequency of multiplicity in X-ray binaries and high-mass 
systems. The current sample size may be insufficient to infer a 
statistically significant relationship, but these results can still 
contribute valuable data to future surveys of high-mass systems. 

This pilot study, in addition to showing evidence for 
potential additional components in X-ray binaries, could 
contribute to stellar multiplicity surveys for massive stars. It 
could also enable us to probe stellar multiplicity at low mass 
ratios (below ~0.1). 


7.2. Stability in Wide Orbit 


The projected separations within the scope of this study 
(from ~350 to ~12,000 au) would suggest that CBCs orbit at 
very large distances from the central X-ray binary. Note that 
CBCs located closer to the X-ray binary within the 2—3 kpc 
distance range cannot currently be detected through direct 
imaging. To assess whether potential CBCs could be 
gravitationally bound to these systems, we calculated the 
binding energy Epina for every source likely to be a CBC. 
Assuming circular orbits, Epina is estimated using the following 
equation (e.g., Naud et al. 2014; Prasow-Emond et al. 2022): 


GMxrpMcomp : (3) 
1.27r 

where G is the gravitational constant, Mxgg is the total mass of 

the central X-ray binary, Mcomp is the mass of the CBC, r is the 


Ebina ~ — 
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projected separation between the CBC and the X-ray binary, 
and 1.27 is the average projection factor between r and the 
semimajor axis assuming a random viewing angle (e.g., 
Brandeker et al. 2006). 

The binding energies for each candidate CBC range from 
roughly —2 x 10? erg to —3 x 10 erg. A binding energy of 
roughly —10"" erg was obtained for the comoving exoplanet 
GU Psc b around an M3 spectral type star (~0.46 Mo), with a 
mass of ^-9Mj,,-13Mjy and located at a distance of ~2000 au 
(Naud et al. 2014). All candidate CBCs have binding energies 
largely exceeding this currently known lower limit. This 
suggests that these sources, if confirmed as CBCs, would fall 
within the gravitational binding range of the X-ray binary. 

We must also consider dynamic stability for systems 
containing more than one candidate CBC. N-body simulations 
have shown that a configuration of two CBCs at the same 
projected separation from the central system can lead to dynamic 
instability (e.g., Kiseleva et al. 1994). Thus, in the case of 4U 
1700—377, where B and C have similar projected separations, 
configurations B+D and C+D are more likely than B+C-+D. 
As for RX J2030.5--4751 and IGR J18483—0311, the candidate 
CBCs have distant projected separations, suggesting that 
configurations B+C and B+C-+D, respectively, are plausible. 


7.3. Companion Formation and Capture Scenarios 


If our findings are confirmed, we hypothesize that CBCs 
orbiting X-ray binaries could originate through two main 
mechanisms: (1) formation within the same environment as the 
central X-ray binary, or (2) capture by the system. On the one 
hand, as detailed in Prasow-Émond et al. (2022), there are three 
scenarios in which CBCs could potentially form in the direct 
environment of X-ray binaries, and these scenarios may unfold 
at different times. First, CBCs could have formed simulta- 
neously with the initial stars that subsequently evolved to form 
the present X-ray binary, resulting from the direct fragmenta- 
tion of a collapsing prestellar core (e.g., Bate 2012). Note that 
this scenario is unlikely because we would expect companions 
with a mass similar to that of the main X-ray binary. Second, 
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Table 4 
Recommended Year for Follow-up Observations for Every X-Ray Binary with 
at Least One Candidate CBC—Calculated Using the Distance, Proper Motions, 
and Astrometric Errors on the 2020 Observations—Enabling a Proper-motion 
Analysis at a >3o Significance Level 


Target Recommended Year 
X Per 2026 
Cyg X-1 2024 
IGR J18483—0311 2024 
1H 22024-501 2027 
4U 1700—377 2024 
RX J2030.5+4751 2024 


their formation might have occurred within the circumbinary 
disk surrounding the initial binary system, prior to the 
supernova explosion of the now-compact object (e.g., Kratter 
et al. 2010). Third, CBCs could have formed within the 
supernova fallback disk arising from the explosion (e.g., 
Wolszczan & Frail 1992). Note that fallback disks do not 
contain enough mass to enable star formation, but they do 
possess enough mass for the formation of substellar objects. 
Additional observations are needed to detect and characterize 
such disks (see Section 7.4). 

On the other hand, stellar and substellar CBCs could be 
gravitationally captured by the X-ray binary. This is analogous 
to the case of the PSR B1620-26 system, where a giant 
exoplanet is thought to have been captured by a binary system 
containing an NS and a WD (e.g., Sigurdsson et al. 2003). 
Given that the HMXBs in our sample are massive (>10 M.) 
and located in the Galactic plane (see Figure 1), such events 
would not be unlikely. 


7.4. Follow-up and Complementary Observations 


The findings of our study primarily consist of intermediate 
results, and additional observations are necessary to conduct 
further analyses and confirm that the CBCs are bound to the 
system. In this section, we provide a list of recommendations 
for follow-up and complementary observations. 

First, we recommend reobserving the systems in the same 
band (L’) using the same instrument (NIRC2) at one or multiple 
additional epochs. This will allow us to conduct proper-motion 
analysis (see Section 6.2) for all the candidate CBCs identified 
in this study. Multiepoch observations will also enable us to 
characterize the orbital motion of these CBCs. Since the 
HMXBs presented in this study are located within 2—3 kpc, the 
time interval between epochs ranges from a few months to a 
few years. Table 4 provides an estimate of the recommended 
year of reobservation for each system, ensuring a statistically 
significant (>30) proper-motion analysis. 

Second, we recommend observing the remaining five 
sources, such as Scorpius X-1, 1A 0620-00, and V404 Cyg, 
to complete the sample of all X-ray binaries within 2—3 kpc 
accessible with Keck/NIRC2. This would also enable us to 
incorporate LMXBs into the analysis and discussion. 

Third, observations in other bands (e.g., K;) would allow us to 
construct color-magnitude diagrams and employ evolutionary 
models to estimate the physical properties of the CBCs with 
greater constraints (e.g., Prasow-Emond et al. 2022). Similarly, 
obtaining the near-infrared spectrum of the candidate CBCs 
would enable us to characterize the nature of the source and 
extract additional physical properties. Finally, submillimeter 
observations of the continuum emission would allow us to detect 
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and characterize potential au-scale circumstellar or protoplane- 
tary disks composed of dust and hot gas (e.g., Coleiro et al. 
2013; Iyer & Paul 2017; Waisberg et al. 2019). 

Contrast curves for all systems with CBCs can be found in 
the Appendix to assess the limit that has been reached during 
these observations. 


8. Summary and Conclusions 


In this study, we presented the first L'-band high-contrast 
images of nearby HMXBs using NIRC2 and the vortex 
coronagraph on the W. M. Keck Observatory. A total of 14 
systems were observed from a sample of 19 X-ray binaries 
within ~2-3 kpc and amenable for direct imaging. One or 
several sources with an S/N > 5 were found in eight of the 
observed X-ray binaries. To discern the nature of these sources 
—whether unrelated objects or candidate CBCs—we employed 
Galactic population models for all systems and proper-motion 
analysis for y Cas. We find that, if confirmed, these results 
would imply the presence of stellar and substellar CBCs in the 
direct environment of X-ray binaries (~350—12,000 au), which 
opens up the discussion on the binary nature of these systems. 
As a pilot study, this work presents a catalog of photometric 
and astrometric parameters for the first epochs of observations. 
Follow-up observations or additional characterization (e.g., 
infrared spectrum) will enable us to conduct proper-motion 
analyses to discriminate more robustly background/foreground 
sources from comoving, gravitationally bound CBCs. 
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Appendix 


This appendix presents additional data related to the main 
text. Figure 6 presents the 5c contrast curves for all X-ray 
binaries with CBCs. Table 5 presents additional relevant 
physical properties for the 14 observed X-ray binaries. Table 6 
presents the properties of the detected sources in the high- 
contrast image, including the status (background/foreground 
source or candidate CBC), optimization parameters, physical 
properties (astrometric and photometric parameters), and 
estimated mass. 


THE ASTROPHYSICAL JOURNAL, 967:8 (19pp), 2024 May 20 


Apparent magnitude (mag) 


Apparent magnitude (mag) 


Apparent magnitude (mag) 


Figure 6. 5c contrast curves in apparent magnitude (left axis) and estimated mass from the age of the X-ray binary (right axis) for y Cas, Cyg X-1, X Per, 4U 1700 
—377, 1H 2202+501, RX J2030.5+4751, and IGR J18483—0311. The curves were generated using VIP for a Gaussian distribution (blue) and with a Student's t-test 


correction (red). 
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Additional Relevant Physical Pd the 14 Observed X-Ray Binaries 

Target Distance PM R.A. PM Decl. m, Age E(B — V) References Lx References Var. Ind. 

(kpc) (mas yr) (mas yr) (+0.5 mag) (Myr) (erg s 1) 
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) 
IGR J18483—0311 2.6 + 0.7 -1.7 40.2 —3.7 € 0.1 7.907 Up to ~30-50 5.22 + 0.02 (3) From ^10? to 10? (3, 4, 5) 
y Cas 0.19 + 0.02 25.7 + 0.5 —3.8 + 0.4 —0.912 8.0 + 0.4 —0.15 (6) ~10°-10°° (7) 
SAX J1818.6—1703 2.3 40.8 —1.6 + 0.2 —4.6 + 0.1 8.964 Up to ~30-50 5.08 + 0.05 (3) From ~10*? to 10% (3, 8) 
1H 2202+501 1.10+0.01 2.36 + 0.01 —0.29 + 0.01 8.285 Up to ~130 0.36 + 0.03 (9) ^9 x 10°? (10) N/A 
4U 22064543 3.1 0.1 —4.17 + 0.02 —3.32 + 0.01 8.700 Up to ~8 0.547 + 0.066 (11) ~10*°-10°6 (12, 13) 
4U 1700-377 1.50.1 2.41 + 0.03 5.02 + 0.02 5.36 ~80 0.50 + 0.01 (14) Up to ~7 x 10? (15) 
IGR J17544—2619 24 02 —0.51 + 0.03 —0.67 + 0.02 7.67 Up to ~8 N/A From ~10* to 10°8 (16, 17) 
RX J2030.5+4751 2.30 + 0.07 —2.71 + 0.02 —4.54 + 0.02 7.088 Up to ~30-50 N/A Up to~10** (18, 19) N/A 
Cyg X-1 2.1 0.1 —3.81 + 0.01 —6.31 + 0.02 6.406 51.5 1.11 + 0.03 (20) ~3 x 107 (21) N 
X Per 0.6+0.1 —1.28 + 0.05 —1.87 + 0.03 4.596 5 ~0.4 (22) ~4 x 10*4 (23) N 
1H 0556+286 150.1 0.63 + 0.03 —2.19 + 0.02 7.618 Up to ~150 N/A Up to ~4 x 10% (24) N/A 
RX J0648.1—4419 0.52 + 0.01 —4.16 + 0.07 5.93 + 0.06 9.150 N/A N/A ~10° (25) N/A 
Vela X-1 2.0+0.1 —4.82 + 0.02 9.28 + 0.02 5.458 Up to ~30-50 ~0.8 (26) Up to ~4 x 10? (27) N 


Note. Column (1): name of the target. Column (2): distance in kpc (see Section 3.3.1). Column (3): proper motion in R.A. in mas yr !. Column (4): proper motion in decl. in mas yr !. Column (5): apparent magnitude in 
the L’ band (see Section 3.3.2). Column (6): estimated age of the X-ray binary in Myr (see Section 3.3.3). Column (7): value of the extinction. Column (8): reference for the extinction. Column (9): observed values of the 
X-ray luminosity in erg s '. Column (10): reference for the X-ray luminosity. Column (11): the variability indicator, as found in Bird et al. (2016). 

References. (1) Schlafly & Finkbeiner 2011; (2) Naze & Motch 2018; (3) Torrejón et al. 2010; (4) Romano et al. 2010; (5) Sguera et al. 2015; (6) Moffat et al. 1973; (7) Raguzova & Popov 2005; (8) Bozzo et al. 2012; 
(9) Simon et al. 2019; (10) Chevalier & Ilovaisky 1998; (11) Nikolov et al. 2017; (12) Negueruela & Reig 2001; (13) Ribó et al. 2006; (14) Hainich et al. 2020; (15) Laurent et al. 1992; (16) Romano et al. 2015; (17) 
Bozzo et al. 2016; (18) Liu et al. 2006; (19) Raguzova 2007; (20) Caballero-Nieves et al. 2009; (21) Di Salvo et al. 2001; (22) Lyubimkov et al. 1997; (23) Coburn et al. 2001; (24) Helfand & Moran 2001; (25) 
Mereghetti et al. 2021; (26) Klare & Neckel 1977; (27) Kreykenbohm et al. 2008. 
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Table 6 


Physical Properties of the Detected Sources in the High-contrast Images 


Target Source Status Ncomp Nov (L! < my) 1 — Pynrelatea(, ©) p 0 my! Est. Mass Proj. Sep. 
(S/N > 5) (%) (mas) (deg) (mag) (au) 
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) 
X Per B cc 20 ~0.1 > 99 587 + 18 112.3 + 2.0 14.4 + 0.6 ~45Mjup—110M yup 350 
Cyg X-1 B cc 26 ~0.6 93 +1 1853 + 12 170.1 + 0.6 16.2 + 0.6 ~0.2—0.3 Mo 4000 
IGR J18483—0311 B cc 2 0.3—0.8 90 +4 2527 + 12 291.5 + 0.2 13.3 + 0.8 ~2.5 Mo 6440 
C cc 3 0.9—1.8 89+4 1516 + 13 21.6 + 0.4 14.8 + 0.8 ~1.1—1.3 Mo 3865 
D bkg 7 3.4—7.0 324 15 2787 + 14 16.9 + 0.2 16.3 + 0.8 es 7100 
E bkg 13 4.1—8.3 14+ 10 3127 + 15 13.8 + 0.2 16.7 + 0.8 7975 
F bkg 17 1.3-3.1 23+ 10 4522 + 14 15.0+0.1 15.3 + 0.8 11530 
G bkg 25 ~7 <5 3761 + 12 69.7 + 0.2 17.0 + 0.8 9600 
H bkg 23 ~ll <5 3358 + 13 83.3 + 0.2 17.6 + 0.8 8560 
I bkg 19 ~ll <5 3168 + 15 106.8 + 0.2 17.6 + 0.8 8080 
SAX J1818.6—1703 B bkg 3 1.5—3.5 25+ 10 2820 + 11 96.9 + 0.2 18.2 + 0.8 6600 
C bkg 20 >30 <5 4069 + 12 353.9 + 0.2 14.9 + 0.8 9525 
1H 2202+501 B cc 14 0.2—0.5 98 +1 1222 + 15 162.4 + 0.5 16.4+0.9 ^ I0My to ~0.4 Mo 1370 
4U 1700—377 B cc 8 0.5—0.7 91 £3 2238 + 12 353.0 + 0.3 14.0 + 0.6 ~1.2-1.3 Mo 4075 
C cc 7 1.3—1.9 68 +4 2616 + 12 274.4 + 0.2 13.4 + 0.6 ~1.5—1.6 Mo 4760 
D cc 3 0.7—0.9 656 4157+ 15 283.4 + 0.1 15.4 + 0.7 7-0.4—0.8 Mo 7570 
IGR J17544—2619 B bkg 5 ~48 Edu) 827 + 12 19.4 + 0.7 16.9 + 0.8 2350 
C bkg 12 ~13 <5 2993 + 13 10.8 + 0.2 15.1 + 0.8 8500 
D bkg 3 ~6 <5 4099 + 11 72.5 + 0.1 13.4+0.8 11640 
E bkg 15 ~86 <5 3170 + 14 70.7 + 0.2 17.3 +0.8 9000 
F bkg 11 ~14 ~T 2392 + 14 103.5 + 0.2 15.3 + 0.8 6790 
G bkg 3 ~22 <5 2947 + 12 118.7 + 0.2 16.1 + 0.8 8370 
H bkg 8 ~8 <5 3935 + 13 147.7+0.1 14.1 + 0.8 11180 
I bkg 10 7-68 «5 2093 + 11 235.1 + 0.3 17.1 + 0.8 5950 
J bkg 15 ~85 <5 3537 + 15 267.1 + 0.2 17.3 + 0.8 10050 
RX J2030.5+4751 B cc 7 ~0.1 > 99 513 +11 302.9 + 0.9 14.8 + 0.8 7-0.3—1.1 Mo 1130 
C cc 17 0.2—0.3 92 £3 3251 + 15 132.1 €: 0.3 16.1 + 0.8 ~0.1—0.6 Mo 7150 
D cc 28 0.3—0.7 7A X8 4208 + 13 343.1 + 0.3 17.1 + 0.8 ~60Mjyp—400M ju, 9250 
y Cas (2017) B 1 2051 + 20 242.8 + 0.3 390 
C 24 1803 + 20 91.4 3: 0.3 343 
^y Cas (2020) B cc 1 2056 + 20 241.2 + 0.3 3.1 + 0.9 ~13 Mo 391 
C bkg 24 1721 + 20 90.8 + 0.3 10.4 + 0.9 ZE 327 


Note. Column (1): the target. Column (2): the label of the detected source (S/N > 5). Column (3): the current status of the source, either background/foreground object (bkg) or candidate CBC (cc). Column (4): the 
optimal number of principal components used for the fitting (comp). Column (5): the expected number of sources in the FOV below the apparent magnitude (npoy (Z^ < mz); see Section 6.1). Column (6): one minus the 
probability of being unrelated with the central X-ray binary (1 — Punrelatea(, ©); see Section 6.1). Column (7): the radial separation from the X-ray binary in mas. Column (8): the position angle in the image in degrees. 
Column (9): the apparent magnitude in the L’ band (no extinction correction was applied). Column (10): the mass of the source, in Mo or Myyp, estimated from evolutionary models (MIST or COND /DUSTY). Column 
(11): the projected separation from the X-ray binary in astronomical units. 
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